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It  is  to  be  noted  here  that  the  first  task  forms  the  major 
portion  of  the  overall  project.  The  work  under  this  task  consisted 
of  developing  computer  program  modules  that  can  be  used  in  general 
purpose  nonlinear  finite  element  programs  to  analyze  the  surface 
instability  of  underground  excavations.  The  objective  of  the 
second,  third  and  fourth  tasks  was  to  assess  the  applicability  of 
our  numerical  approach  in  the  case  of  a  simple  boundary  value 
problem.  The  numerical  analysis  had  to  be  calibrated  for  a  simple 
boundary  value  problem  before  extending  it  to  the  complex 
geometries  and  material  properties  encountered  in  the  case  of  real 
excavations.  The  calibration  test  chosen  is  a  plane  strain 
compression  of  rock  sample.  It  will  be  referred  to  as  the  wedge 
test.  During  the  second  task,  the  program  modules  developed  in  the 
first  task  are  used  to  obtain  numerical  results  in  the  particular 
case  of  the  wedge  test.  The  performance  of  these  finite  element 
modules  are  tested  by  comparing  their  results  with  analytical 
results  or  other  simpler  numerical  solutions.  These  additional 
results  were  obtained  by  using  the  program  BIF,  which  was 
especially  written  to  carry  out  the  parametric  study  under  the 
third  task.  The  parametric  study  consisted  of  investigating  the 
influence  of  the  parameters  of  various  constitutive  equations  on 
the  buckling  load  of  the  wedge  test.  The  fourth  and  last  task  was 
to  investigate  how  much  surface  instability  was  influenced  by  the 
introduction  of  friction  forces  at  the  boundaries  of  the  wedge 
test. 

The  numerical  method  which  has  been  developed  in  this  report 
has  only  been  applied  to  analyze  surface  instability  of  the  wedge 
test  but  not  of  real  excavations.  The  extension  of  the  method  to 
real  excavations,  which  are  more  complex  boundary  value  problems, 
is  the  long  term  objective  of  the  project.  This  extension  will  be 
the  object  of  future  research. 
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1.2  REPORT  ORGANIZATION. 


The  main  body  of  this  report  is  comprised  of  5  major  sections 
covering  respectively 

1.  the  definition  and  past  work  on  rockbursting  and  surface 
instability, 

2.  the  finite  element  techniques  used  to  describe  the  phenomenon 
of  surface  instability, 

3.  the  numerical  results  obtained  in  the  particular  case  of  the 
wedge  test, 

4.  the  influence  of  material  modeling  on  surface  instability,  and 

5.  the  influence  of  friction  boundary  conditions  on  surface 

instability. 

Mathematical  derivations  have  been  documented  in  appendices. 

1.3  DEFINITION  OF  ROCKBURSTING. 

The  term  rockburst  is  used  to  designate  a  violent  failure  of 
rock  which  is  sometime  experienced  in  deep  underground 
excavations.  It  involves  a  rapid  convergence  and  oscillation  of 
the  excavation  walls,  followed  by  slabbing  and  failure  of  the  rock 
immediately  adjacent  to  the  excavation  (Ortlepp,  1978) .  This 
violent  release  of  energy  has  been  observed  on  a  scale  ranging 
from  the  expulsion  of  small  fragment  of  rock  to  major  collapse  of 

the  excavation.  It  is  often  observed  to  follow  enlargement  of  the 

cavity  by  blasting.  Rockburst  appears  to  be  more  frequent  in  rock 
which  are  hard  and  brittle. 

The  nature  and  the  circumstances  accompanying  this  phenomenon 
suggest  that  underground  excavations  in  hard  rock  subject  to 
dynamic  loading  caused  by  nuclear  explosion  are  prone  to 
rockbursting.  In  mining  operations  the  mitigation  of  rock  burst 
relies  on  designing  the  layout  of  the  mine  excavation  and  the 


sequence  of  extraction  in  order  to  reduce  as  much  as  possible  the 
vertical  stress  concentration  around  the  mine  (method  of  energy 
release  rate,  Jaeger  and  Cook,  1977,  Salamon,  1974).  Although  the 
method  provides  a  practical  design  tool  for  reducing  the  hazards 
of  rockburst,  it  does  not  provide  a  mean  of  calculating  whether 
rockburst  will  occur  or  not. 

1.4  PAST  WORK  ON  ROCKBURSTING. 


m 
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Several  phenomenological  models  have  been  proposed  for 
explaining  rockbursting  phenomena.  They  are  based  upon  limiting 
static  equilibrium,  nonlinear  elasticity,  strain-softening 
material,  unstable  propagation  of  pre-existing  cracks,  and  finally 
surface  instability.  Lippman,  1978,  Avershin  et  al.,  1972,  examine 
the  static  equilibrium  of  rock  block  adjacent  to  the  excavation. 
The  rockburst  occurs  in  circumstances  where  the  resisting  forces 
do  not  balance  the  lateral  forces  tending  to  propel  the  block 
inside  the  cavity.  Freudenthal,  1977,  use  nonlinear  elasticity 
(hyperelasticity)  to  account  for  shear  induced  volume  expansion. 
He  considered  that  rockburst  and  spalling  of  the  opening  walls  are 
caused  by  crack  propagation  associated  with  the  tensile  stress 
around  the  excavation.  Salamon,  1974,  and  Pietruszczak  and  Mroz, 
1980,  Pariseau,  1979,  characterize  the  rock  as  a  strain-softening 
material.  Instability  becomes  possible  once  the  rock  has 
experienced  its  peak  strength.  This  mechanism  can  however  be 
refuted  on  the  ground  that  the  softening  observed  in  compression 
of  brittle  rock  specimen  is  test  dependent  and  is  always 
associated  with  localized  failure  mechanisms  (Drescher  and 
Vardoulakis,  1982).  Nemat-Nasser  and  Horii,  1982,  analyze  rock 
bursting  as  buckling  of  slabs.  The  slabs  become  detached  from  the 
rock  by  tensile  extension,  parallel  to  the  free  surface  of 
pre-existing  cracks,  which  are  caused  by  far-field  compressive 
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stress.  Finally,  Vardoulakis,  1984,  uses  the  theory  of  bifurcation 
to  analyze  rockburst  as  surface  instability. 

1.5  PAST  WORK  ON  SURFACE  INSTABILITY. 

The  pioneering  work  on  surface  instability  was  carried  out  by 
Biot,  1965,  who  examined  the  half-plane  problem  with 
incompressible  isotropic  and  anisotropic  elastic  material.  A 
simple  example  of  the  surface  instability  treated  by  Biot  in  1965 
is  summarized  in  Appendix  A.  In  this  example,  Biot  showed  how  the 
elastic  half-space  looses  its  rigidity  for  some  value  of 
compressive  stress  and  that  this  loss  of  rigidity  corresponds 
to  the  phenomenon  of  surface  buckling  or  surface  instability.  The 
buckling  modes  of  the  half -space  are  surface  waves,  the  amplitude 
of  which  decays  rapidly  with  depth.  Hutchinson  and  Tvergaard, 
1980,  extended  Biot's  results  for  elastoplastic  materials.  They 
noted  that  the  existence  of  surface  buckling  strongly  depends  on 
the  type  of  plasticity  which  is  employed  to  describe  material 
behavior.  The  instability  was  found  to  occur  at  low  loads  for 
solids  characterized  by  a  deformation  theory  of  plasticity  and  at 
high  and  unrealistic  loads  for  solids  characterized  by  the  flow 
theory  of  plasticity.  Vardoulakis,  1984,  basing  his  analysis  on  a 
deformation  of  plasticity  showed  that  the  bifurcation  load  marks 
the  beginning  of  dynamic  instable  conditions  and  that  the 
existence  of  surface  buckling  is  related  to  the  ratio  of  the 
tensile  over  the  compressive  strength.  This  ratio  can  be  used  as  a 
measure  of  the  brittleness  of  the  rock;  the  smaller  the  ratio,  the 
more  brittle  the  rock  is.  Horii  and  Nemat-Nasser ,  1982,  departing 
from  the  assumption  of  incompressibility  made  in  the  previous 
investigations,  introduced  dilatancy  and  showed  that  the  inclusion 
of  this  material  property  reduces  the  magnitude  of  the  bifurcation 
load.  Horii  and  Nemat-Nasser  also  noted  that  the  occurrence  of 
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NUMERICAL  METHOD  FOR  ASSESSING  SURFACE  INSTABILITY 

The  following  section  summarizes  the  numerical  methods  which 
have  been  used  in  the  numerical  assessment  of  rockbursting.  This 
section  only  outlines  major  results.  Mathematical  derivations  are 
documented  either  in  appendices  or  in  the  publications  referenced 
in  the  text.  The  methods  presented  in  this  section  have  been 
implemented  in  the  nonlinear  finite  element  code,  LINOS.  The 
computer  code  LINOS  was  developed  at  the  University  of  Southern 
California  three  years  ago.  Its  main  objective  is  to  serve  as  a 
vehicle  for  applying  nonlinear  constitutive  equations  of  soils  and 
rocks  especially  in  the  range  of  large  deformation  and  in  presence 
of  bifurcations.  LINOS  was  therefore  especially  suited  for  the 
analysis  of  surface  instability.  According  to  the  modular 
structure  of  LINOS,  which  is  common  to  most  of  the  large  finite 
element  codes,  the  library  of  element  and  material  subroutines  can 
be  easily  expanded  without  alteration  of  the  main  code.  In  the 
context  of  this  project,  a  new  element  module  and  two  new  material 
modules  were  added  to  the  element  and  material  library.  In 
addition  to  these  three  new  modules,  two  other  modules  have  been 
coded  to  perform  the  calculation  of  the  eigenvalues  and 
eigenvectors  of  the  stiffness  matrix.  All  these  modules  can 
certainly  be  modified  to  fit  into  other  finite  element  codes  such 
as  ADINA  (Bathe,  1982) . 

This  section  is  divided  into  six  subsections: 

1.  updated  Lagrangian  variational  formulation 

2.  calculation  of  tangential  stiffness  matrix 

3.  detection  of  bifurcation  point 


4.  calculation  of  eigenvalues  and  eigenvectors  of  tangential 

stiffness  matrix 

5.  constitutive  models 

6.  finite  element  implementation  of  constitutive  models 

2.1  UPDATED  LAGRANGIAN  VARIATIONAL  FORMULATION. 

Biot,  1965,  Horii  and  Nemat-Nasser ,  1982,  and  Vardoulakis, 
1984,  based  their  analysis  of  surface  instability  upon  the  partial 
differential  equations  of  equilibrium  for  incremental  stresses. 
They  did  not  use  the  variational  principles  which  are  the 
fundamental  basis  of  finite  element  methods.  One  of  the  first  and 
most  important  point  of  the  numerical  assessment  of  surface 
instability  was  to  verify  the  compatibility  of  the  nonlinear 
variational  principle  of  our  finite  element  code,  LINOS,  with  the 
nonlinear  partial  differential  equations  used  by  previous 
investigators.  This  compatibility,  which  is  rather  straightforward 
to  derive  in  the  case  of  small  deformation,  becomes  less  obvious 
in  the  case  of  large  deformation  due  to  the  multiple  choices  of 
accounting  for  material  and  geometrical  nonlinearities.  In  fact 
the  variational  principle  used  in  LINOS  is  identical  to  the 
variational  principle  derived  by  Biot  (1965,  pp  73-79),  and  by 
McMeeking  and  Rice,  1978.  This  formulation,  which  requires  the 
update  of  the  body  configuration  as  it  deforms,  is  referred  to  as 
updated  Lagrangian.  The  position  of  the  finite  element  mesh 
changes  after  each  time  step  in  order  to  perform  the  volume  and 
surface  integrations  on  the  deformed,  not  initial,  body 
configuration.  The  updated  Lagrangian  method  is  detailed  by  Bathe, 
1982,  (pp  335-406) . 


2.2  CALCULATION  OF  TANGENTIAL  STIFFNESS  MATRIX. 


The  finite  element  analysis  of  bifurcation  problems  requires 
more  accurate  and  numerous  calculation  of  the  tangential  stiffness 
matrix  than  most  conventional  nonlinear  finite  element  analysis. 
In  nonlinear  computations,  accurate  results  can  be  obtained  with 
an  approximate  stiffness  matrix  if  an  iterative  technique  balances 
the  internal  and  external  forces.  For  example,  in  the  case  of  the 
modified  Newton-Raphson  method,  the  tangential  stiffness  matrix 
can  be  reformed  only  once,  at  the  first  time  step,  or  a  few  times, 
at  some  later  time  steps.  For  modified  Newton-Raphson,  iterations 
are  carried  out  until  the  out-of-balance  or  residual  forces  become 
negligible.  Out-of-balance  or  residual  forced  are  equal  to  the 
difference  between  internal  and  external  forces.  Modified 
Newton-Raphson  methods  are  sometimes  preferred  to  full 
Newton-Raphson  technique  since  they  are  less  time-consuming,  i.e., 
the  stiffness  matrix  does  not  have  to  be  reformed  as  many  times  as 
for  full  Newton-Raphson. 

However,  in  the  case  of  bifurcation  analysis,  the  stiffness 
matrix  must  be  calculated  frequently  and  accurately.  The  objective 
of  a  finite  element  analysis  of  bifurcation  problems  is  to  detect 
the  loss  of  uniqueness  of  the  incremental  solution  Ad  for  the 
following  general  type  of  nonlinear  matrix  problem  of  n  dimension 

K(d) .Ad  =  AF  (2) 

where 

Ad  =  increment  of  nodal  displacement  at  time  t 

d  =  nodal  displacement  at  time  t 

(cumulated  value  of  Ad  between  times  0  and  t) 
aF  =  residual  or  out-of-balance  force 

(difference  between  internal  and  externally  applied  loads) 
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K(d)=  tangential  stiffness  matrix 

(dependent  on  the  nodal  displacement  d  at  time  t) 
n  =  number  of  degrees  of  freedom 

The  bifurcation  of  the  increment  Ad  from  one  to  another  type 
of  solution  is  only  possible  when  the  problem  (2)  admits  more  than 
one  solution.  This  loss  of  uniqueness  can  only  take  place  when  the 
matrix  K(d)  becomes  singular,  which  is  associated  with  a  zero 
determinant.  Therefore,  the  detection  of  possible  points  of 
bifurcation  requires  an  accurate  evaluation  of  the  stiffness 
matrix,  especially  in  the  vicinity  of  the  bifurcation.  In  order  to 
compute  the  stiffness  matrix  accurately,  a  new  element  has  been 
developed.  This  element  is  a  two-dimensional  plane-strain 
isoparametric  four-node  element  with  four  stress  point  which  can 
handle  large  deformation.  The  four  stress  points  are  located  at 
the  gaussian  points  of  the  element.  One  stress  point  can  also  be 
selected  at  the  element  centroid.  The  integration  rules  which  are 
used  to  evaluate  the  volume  integrals  at  the  element  level  can  be 
chosen  with  two  by  two  or  one  Gaussian  points.  In  the  present 
work,  only  four  stress  points  and  four  integration  points  were 
used.  The  influence  of  reduced  integration  and  stress  points  on 
the  surface  instability  has  not  been  investigated  and  could  be  the 
object  of  future  research. 

2.3  DETECTION  OF  BIFURCATION  POINT. 
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The  bifurcation  point,  which  corresponds  to  the  loss  of 
uniqueness  of  the  incremental  matrix  problem  (2) ,  is  possible  when 
the  determinant  of  the  increment  tangential  matrix  becomes  equal 
to  zero.  In  the  present  numerical  approach,  the  bifurcation  point 
can  be  detected  by  choosing  one  out  of  three  techniques.  It  can  be 
found  by  checking  the  sign  change  of  either  the  pivots,  the 
determinant,  or  the  eigenvalues  of  the  stiffness  matrix.  The  first 
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method  of  detection  is  rather  simple  and  efficient.  It  takes  place 
right  after  the  factorization  of  the  stiffness  matrix  into  upper, 
lower  triangular  and  diagonal  matrices: 


K(d)  =  L.D.U 


L  = 


D  = 


U  = 


lower  triangular  matrix 
diagonal  matrix 
upper  triangular  matrix 


By  definition,  the  pivots  of  K  are  the  diagonal  entries  of  D, 
noted  i=l,n.  In  the  case  of  symmetric  matrix, 


L  =  U 


where  U  represents  the  transpose  of  U. 

The  second  method  computes  the  determinant  of  the  stiffness  matrix 
K  by  calculating  the  products  of  the  pivots.  Since  for  triangular 


matrices  L  and  D 


det(L)  =  det(U)  =1 


then  the  determinant  of  K  is  the  products  of  the  pivots; 


det(K)  =  det(D)  =  .  0^2  •  .  . 


The  third  and  last  method  is  more  time  consuming  than  the  first 
two  methods;  it  requires  the  calculation  of  the  eigenvalues  of  the 
stiffness  matrix.  This  calculation  is  carried  out  just  after  the 
formation  of  the  tangential  stiffness  matrix  but  before  its 


factorization.  If  the  selected  bifurcation  condition  is  met  at 


time  t„ . ,  ,  then  bifurcation  takes  place  between  times  t  and  t 

n+i'  ^  n  n+1 


**  -* 
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If  the  bifurcation  time  is  needed  with  more  accuracy,  the 
step-by-step  analysis  can  be  restarted  from  time  t^^  by  using  a 
smaller  time  interval,  and  stopped  when  the  bifurcation  condition 
is  met  again.  Theoretically,  this  process  could  be  repeated 
several  times  until  the  time  interval  [t  ,t  becomes  small 
enough.  Alternate  techniques  of  detection  use  Regula-Falsi  method 
or  the  method  proposed  by  M.  Fu j ikake  (1985).  Iterative  techniques 
to  refine  the  bifurcation  time  was  not  required  in  the  present 
analysis  since  very  small  loading  steps  were  used. 

2.4  CALCULATION  OF  EIGENVALUES  AND  EIGENVECTORS  OF  THE  TANGENTIAL 
STIFFNESS  MATRIX. 

Two  program  modules  have  been  coded  to  compute  the 
eigenvalues  and  eigenvectors  of  the  tangential  stiffness  matrix. 
One  modules  calculates  only  the  eigenvalues  during  the 
step-by-step  integration  in  order  to  detect  the  bifurcation  point. 
The  second  module  calculates  both  the  eigenvalues  and  the 
eigenvectors  of  the  stiffness  matrix  after  the  detection  of  the 
bifurcation  point.  With  the  help  of  the  interactive  capability  of 
the  second  module,  the  user  can  select  and  plot  the  eigenvectors 
as  it  is  important  to  understand  the  significance  of  the 
eigenvectors  corresponding  to  the  lowest  eigenvalues.  Both  modules 
use  subroutines  from  the  EISPACK  package  (Smith,  1976) . 

2.5  CONSTITUTIVE  MODELS. 


In  the  context  of  this  analysis,  the  following  general  class 
of  rate-type  constitutive  equations  are  considered: 


where  r  is  the  Jauman  rate  of  Kirchhoff  stress,  D  the  rate  of 
deformation,  and  C.  ,  may  be  any  function  of  states  of  stress 

X  ^  iC  X 

and  strain  which  characterizes  the  incremental  material  response. 
In  the  following  section,  three  types  of  constitutive  equations 
(7)  will  be  specified:  isotropic  and  anisotropic  hypo-elastic 
materials,  and  isotropic  elasto-plastic  materials  of  the  flow 
theory  of  plasticity.  The  selected  constitutive  models  have  two, 
three,  or  five  material  constants,  respectively,  which  quantify 
the  response  of  a  particular  material. 

2.5.1  Isotropic  hypo-elastic  solids:  (two  material  constants). 
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.  i  ,  j  -1 , 2 , 3 

sum  on  k, 1=1, 2, 3 


(8) 


where  S is  Kroneker  symbol; 

6ij  =  1  if  i=j,  and  6^^  =0  if  i^j  i,j=l,2,3  (9) 

and  G  and  \  are  shear  and  Lame's  moduli,  respectively. 

2.5.2  Anisotropic  hypo-elastic  solids  (three  material  constants). 


As  a  particular  case  of  anisotropy,  transverse  anisotropy  in 
plane  strain  can  be  accounted  for  by  introducing  a  transverse 
shear  modulus  G^.  In  the  case  of  incompressible  material  with  the 
plane  of  isotropy  normal  to  the  direction,  the  plane  strain 
constitutive  equation  is; 
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incompressibility 
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deviatoric  stress  rate 
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deviatoric  strain  rate 

(lOd) 

A  new  material  model  is  introduced  by  modifying  relation 
order  to  have  compressibility; 

(10)  in 
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This  material  model,  which  is  not  transversely  isotropic  in  a 
rigorous  sense,  is  referred  to  as  anisotropic  hypo-elastic.  If 
the  shear  moduli  are  equal 

G  =  (11c) 
then  the  equation  (11)  is  the  one  of  the  fully  isotropic  case. 

2.5.3  Elastoplastic  solid  (five  material  constants). 

Adopting  the  notation  of  Rudnicki  and  Rice  (1975),  the 
following  elastoplastic  constitutive  relations  of  the  flow  theory 
of  plasticity  are  considered: 
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where 


G  =  shear  modulus 

B  =  bulk  modulus 

H  =  plastic  modulus 

M  =  friction  coefficient 

0  =  dilatancy  parameter 
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deviator  stress 
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invariant  of  deviator  stress  (12c) 
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mean  pressure 


(12d) 


As  shown  in  the  stress  space  (a,T)  of  Fig.l,  the  friction 
coefficient  n  and  dilatancy  parameter  p  are  the  slope  of  yield 
surface  and  normal  of  plastic  potential  surface,  respectively. 
The  relation  (12. a)  is  isotropic  and  depends  only  on  two  stress 
invariants,  a  and  r.  They  are  independent  of  the  third  stress 
invariant,  which  is  sometime  referred  to  as  Lode's  angle. 
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Figure  1.  Geometrical  interpretation  of  the  coefficient 
of  friction  y  and  the  dilatancy  factor  g 
(after  Rudnicki  and  Rice,  1978). 

The  deformation  theory  of  plasticity  has  not  been  considered 
in  this  analysis,  although  it  was  found  to  yield  more  realistic 
prediction  of  buckling  load  than  the  flow  theory  of  plasticity 
(Hutchinson  and  Tvergaard,  1980,  Vardoulakis,  1984) .  The 
superiority  of  the  deformation  theory  over  the  flow  theory  in  the 
particular  case  of  buckling  analysis  results  from  the  dependence 
of  the  incremental  response  upon  the  stress  rate  direction. 
However,  in  some  models  of  the  deformation  theory,  the 
constitutive  matrix  is  a  discontinuous  function  of  the  stress 
increment  direction.  Such  a  discontinuity  may  cause  numerical 
problems  in  the  context  of  finite  element  analyses.  Therefore, 
the  deformation  theory  was  not  considered  for  this  reason. 

A  particular  case  of  the  general  elastoplastic  relation  has 


o 


been  specified  by  Horii  and  Nemat-Nasser  (1982) .  Their  model 
describes  the  particular  experimental  data  obtained  on 
sandstones.  The  plastic  modulus  H  is  a  function  of  the  first 
stress  invariant  a  and  the  effective  strain  7.  The  coefficient  ^ 
is  a  constant,  while  the  dilatancy  angle  $  depend  on  the  stress 
states  (see  Horii  and  Nemat-Nasser,  1982,  for  further  detail) 


2.6  FINITE  ELEMENT  IMPLEMENTATION  OF  CONSTITUTIVE  MODEL. 


From  the  point  of  view  of  nonlinear  finite  element  methods 
the  constitutive  models  are  required  in  two  operations: 

(1)  formation  of  the  stiffness  matrix  and  (2)  calculation  of  the 
residual  forces.  Consequently  the  material  modules  must  perform 
two  main  tasks:  (1)  calculation  of  the  constitutive  equation 
and  (2)  stress  calculation.  The  first  task  is  done  by 
evaluating  the  coefficients  as  specified  by  the  relations 

(8),  (9)  and  (10).  The  second  task  is  less  straightforward  than 

the  first  task.  The  strain  increment  a«  and  rigid  rotation  R 
which  result  from  incremental  nodal  displacement  are  not 
infinitesimal  but  finite.  Appropriate  numerical  techniques  are 
therefore  required  to  integrate  the  rate  equations  (8) ,  (9)  and 
(10)  between  time  t^  and  This  integration  is  performed  by  a 

module  which  is  referred  to  as  a  stress  calculation  subroutine. 
Stress  calculation  subroutines  play  an  important  role  in  the 
convergence  of  the  iterative  scheme  which  balances  out  internal 


and  external  forces.  If  the  stress  are  not  calculated  with 


sufficient  accuracy,  convergence  may  be  hard  to  achieve.  This 
absence  of  convergence  is  especially  to  be  expected  in  the 
vicinity  of  bifurcation  point,  where  the  stiffness  matrix  is 
about  to  become  singular.  Error  in  the  stress  calculation 
subroutine  may  generate  artificial  out-of-balance  force  which 
prevents  the  system  from  reaching  a  bifurcation  point. 


I 
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In  the  case  of  finite  deformation  the  stress  increment 
calculation  are  carried  out  by  using  Hughes  and  Winget  (1981) 
method.  The  objective  stress  rate  of  relation  (7)  involves  the 
spin  W  in  order  to  account  for  the  rigid  body  rotations  of  a 
material  point  relative  to  the  spatial  coordinates.  In 
incremental  form,  instead  of  rate  form,  the  spin  W  becomes  an 
orthogonal  rotation  through  an  incremental  angle.  Hughes  and 
Winget  (1981)  have  given  a  modern  account  of  this  process  and 
have  provided  a  direct  way  to  evaluate  the  orthogonal  rotation 
matrix  R  from  the  spin  W.  Thus 

R  =  (5  -  AtW/2)~^.(5  +  AtW/2)  (13 


where  S  is  Kroneker's  symbol,  and  At  the  time  interval.  In 
two-dimension  R  can  be  written  as: 


cose 


-  sintf 


sin^ 


cose 


.  .  1/2 
Half  angle  trigonometric  formulas  are  used  to  get  R  '  the  square 

root  of  R.  With  these  constructions  the  rate-type  constitutive 

models  (7)  can  be  integrated  over  the  increment  from  t  to  t  , , 

n  n+l 

in  the  following  way.  First  the  stresses  which  represent  the 
stress  at  all  the  stress  points  at  time  t^,  are  advanced  to  the 

Vi/2  =  (Wi)/2  by 

n+1/2  ^  n' ^ 

Using  the  intermediate  configuration  between  time  t^  and 

and  the  time  step  At,  the  strain  increments  ^*n+l/2 

calculated  at  all  the  stress  points  from  the  nodal  displacement 

between  times  t^  and  Then  the  constitutive  equation  is 

integrated  and  new  stresses  a  are  obtained.  These  are  then 

n+l/2 


rotated  from  t  , ,  to  t  . ,  by  the  same  process  as  in  equation 

n+1/2  n+l  •* 

(14).  The  integration  of  the  constitutive  equation  is  carried  out 

by  using  a  subincremental  technique  dividing  the  prescribed 

strain  increment  Ae  ,  ,  into  several  subincrements.  Each 

n+l/2 

subincrement  is  integrated  with  an  improved  Euler  method  without 
iteration.  In  our  computations,  the  subincrementation  was  not 
fully  exploited;  only  one  increment  was  used. 


SECTION  3 


NUMERICAL  RESULTS  ON  SURFACE  INSTABILITY  OF  THE  WEDGE  TEST 


The  following  section  summarizes  the  numerical  results  which 
are  obtained  on  the  surface  instability  of  a  particular  boundary 
value  problem  -  the  plane-strain  compression  of  a  cube,  which  is 
referred  to  as  the  wedge  test.  The  calculations  were  carried  out 
by  using  computer  modules  which  are  based  on  the  methods  and 
principles  presented  in  the  previous  sections.  The  objectives  of 
the  present  section  are  (1)  to  verify  that  our  computer  modules 
are  correctly  implemented  and,  (2)  to  check  that  the  selected 
finite  element  method  is  appropriate  to  describe  surface 
instability  in  the  particular  case  of  the  wedge  test.  The 
performance  of  the  computer  modules,  especially  the  new  element 
and  material  model  subroutines,  are  tested  by  comparing  finite 
element  results  with  analytical  or  numerical  results  derived  for 
hypo-elastic  and  elastoplastic  material  models. 

This  section  is  divided  into  six  subsections: 

1.  Geometry  and  boundary  conditions  of  the  wedge  test 

2.  Finite  element  mesh 

3.  Analytical  solutions  for  stress-strain  response  and 
bifurcation  load 

4.  Integration  algorithm 

5.  Change  of  eigenvalues  versus  applied  displacement 

6.  Eigenvalues  and  eigenvectors  at  the  bifurcation  point 

3.1  GEOMETRY  AND  BOUNDARY  CONDITIONS  OF  THE  WEDGE  TEST. 

The  geometry  and  boundary  conditions  of  the  wedge  test  are 
shown  in  Fig. 2.  The  cubical  body  has  a  one  to  one  aspect  ratio. 


The  initial  length  of  its  sides  is  equal  to  unit  length.  The  left, 
right,  and  lower  boundaries  are  frictionless.  The  upper  boundary 
is  free  of  traction.  The  block  is  subjected  to  plane  strain 
deformation  by  preventing  displacement  in  the  direction.  The 
horizontal  displacement  of  the  nodes  located  on  the  left  boundary 
are  prescribed  while  their  vertical  displacement  is  free.  The 
horizontal  displacement  u^  on  the  left  boundary  is  specified  by 
applying  a  selected  displacement  increment  Au^,  the  size  of  which 
is  0.0025.  Therefore  200  increments  are  required  to  reach  the 
horizontal  displacement  of  0.5.  In  this  deformation  range  -  50%  - 
the  strain  is  finite,  which  results  in  nonlinear  relationships 
between  strain  and  displacement  gradient.  Such  a  geometrical 
nonlinearity  constitutes  a  meaningful  test  for  the  updated 
Lagrangian  technique.  The  vertical  nodal  displacements  on  the  left 
boundary  will  be  fixed  in  a  later  section  to  investigate  the 
effect  of  nonuniform  stress  field  on  surface  instability.  The 
geometry  and  boundary  conditions  of  Fig. 2  correspond  to  the 
compression  of  a  material  block  between  five  lubricated  rigid 
plates  with  only  one  free  surface.  For  such  an  idealized  test 
condition,  the  trivial  solution  is  a  uniform  state  of  stress  and 
strain  within  the  sample  even  for  large  deformation. 

3.2  FINITE  ELEMENT  MESH. 


Four  different  meshes,  with  one,  four,  nine  and  twenty  five 
elements  respectively,  are  used.  Fig. 2  shows  only  the  25  elements 
model.  All  the  elements  for  a  given  mesh  have  the  same  initial 
size.  They  are  four  nodes  isoparametric  elements  with  four  stress 
points  per  element  located  at  the  Gaussian  points  of  the  element. 
The  number  of  degrees  of  freedom  of  the  problem  varies  with  the 
number  of  elements  as  indicated  in  Table  1.  As  shown  in  Table  1, 
the  number  of  nodal  degrees  of  freedom  increases  rapidly  with  the 


number  of  elements,  which  causes  a  rapid  increase  of  the 
computation  time  to  extract  the  eigenvalues  and  eigenvectors, 
especially  in  the  case  of  matrices  which  are  neither  symmetric  nor 

positive  definite.  Mainly  for  this  reason,  the  number  of  elements 

2 

n  was  limited  to  25.  Higher  order  elements  or  nonuniform  mesh  are 
not  used  in  this  project.  A  more  detailed  analysis  of  the  effect 
of  the  mesh  refinement  on  surface  instability  need  to  be  carried 
out  in  the  future.  Surface  instability,  which  Biot  associates  to 
the  emergence  of  surface  waves  the  amplitude  of  which  decays 
exponentially  with  depth,  could  certainly  be  described  more 
accurately  by  refining  the  mesh  in  the  vicinity  of  the  free 
surface . 


Table  1.  Degree  of  freedom  versus  number  of  elements  in  the  wedge 
test 


number  of  elements 
(on  1  side)  (total) 


number  of  degrees  of  freedom 
(size  of  stiffness  matrix) 


3.3  ANALYTICAL  SOLUTION  FOR  STRESS-STRAIN  RESPONSE  AND  BIFURCATION 
LOAD. 


The  idealized  hypo-elastic  model  with  constant  moduli  was 
considered  in  order  to  obtain  simple  analytical  solutions.  This 
material  may  not  represent  the  behavior  of  any  real  material  but 
it  leads  to  simple  analytical  solutions,  even  for  large 
displacement,  when  compared  to  other  realistic  models.  The 
analytical  derivation  is  detailed  in  the  Appendix  B.  The  material 


parameters  of  hypo-elastic  models  are  listed  in  Table  2.  The 
buckling  stress  of  the  free  surface  of  the  half-space  made  of  such 
an  idealized  material  can  also  be  calculated  analytically  or 
numerically.  The  values  of  stress  and  displacement  at  the  onset  of 
surface  instability  are  shown  in  Table  3  and  Fig. 3.  The 
stress-displacement  response  of  Fig. 3  calculated  by  the  finite 
element  analysis  of  hypo-elastic  material  agree  perfectly  with  the 
analytical  results  described  in  the  appendix  B,  which  demonstrates 
that  the  element  and  material  model  were  correctly  implemented. 

Table  2.  Material  constants  of  hypo-elastic  and  elastoplastic 
models 

model  material  constant 

isotropic  hypo-elastic  shear  modulus  G  1. 

Poisson  ratio  0.3 

anisotropic  hypo-elastic  shear  modulus  G  1. 

transverse  modulus  G.  0.1 

Poisson  ratio  w  0.3 

elastoplastic  Young's  modulus  E  200. 

Poisson  ratio  i/  0.3 

Friction  coefficient  ^  0.39 

variable  dilatancy  0 

0  =  J3n  +  r/a 

variable  plastic  modulus  H 
H  =  exp(1.8-0.36a+e°'®~^°^-227) 


Table  3.  Analytical  values  of  displacement  and  normalized 
stress  at  the  buckling  of  the  half  space  for 

hypo-elastic  and  elastoplastic  models. 


material  model 

^1 

^,l/G 

isotropic 

-0.3593 

-1.4484 

anisotropic 

-0.07336 

-0.2225 

elastoplastic 

-0.058 

-0 . 0334 

It  is  important  to  notice  that,  for  the  particular  case  of  the 
anisotropic  hypo-elastic  material  subjected  to  the  wedge  test,  the 
transverse  shear  modulus  affects  only  the  buckling  stress  but 
not  the  stress-strain  response.  Isotropic  and  anisotropic 
hypo-elastic  material  give  the  same  stress-strain  response  but 
different  buckling  stress.  Therefore  a  good  fit  of  the 
stress-strain  response  of  the  wedge  test  does  not  necessarily  mean 
a  correct  prediction  of  the  buckling  load  since  surface  buckling 
depends  on  an  aspect  of  the  material  behavior  which  is  not 
apparent  in  the  stress-strain  response. 

In  the  case  of  the  elastoplastic  model  of  Horii  and 
Nemat-Nasser ,  1982,  the  stress-displacement  of  Fig. 4  has  been 
obtained  by  integrating  step-by-step  the  plane  strain  constitutive 
equation  for  the  material  constants  specified  in  Table  2  (see 
Horii  and  Nemat-Nasser,  1982,  for  detail).  Surface  instability  was 
detected  by  using  a  technique  similar  to  the  one  employed  in  the 
computer  code  BIF.  The  values  of  displacement  and  normalized 
stress  are  shown  in  Table  3. 

It  is  worth  noting  that,  for  both  hypo-elastic  and 


the  block  of  finite  size  shown  in  Fig. 2,  and  that  the  buckling 
stress  was  obtained  for  the  infinite  half-space. 

3.4  INTEGRATION  ALGORITHM. 

Two  integration  schemes  were  alternately  used  to  integrate 
step-by-step  the  nonlinear  finite  element  problem  of  the  wedge 
test;  full  and  modified  Newton-Raphson  methods.  During  the 
modified  Newton-Raphson  computations,  the  stiffness  matrix  was 
reformed  only  once  each  time  step,  which  was  less  time  consuming 
than  full  Newton-Raphson,  but  accurate  enough  due  to  the  small 
size  of  the  loading  increment.  At  a  prescribed  boundary 
displacement  u^,  the  convergence  of  both  methods  is  measured  by 
the  norm  of  the  residual  load  vector  after  a  certain  number  of 
iterations.  Iterative  calculations  stop  when  the  norm  of  the 
residual  force  is  less  than  some  given  convergence  parameter  TOL. 
Generally  TOL  is  about  equal  to  lo”^.  Such  a  value  was 
sufficiently  small  to  maintain  a  uniform  stress  and  strain  state 
within  the  sample  of  the  wedge  test,  and  to  obtain  the  bifurcation 
point.  In  some  circumstances,  larger  values  of  TOL  were  found  to 
lead  to  heterogeneous  solutions.  The  error  in  the  computed 
solution  has  no  important  effect  far  away  from  the  bifurcation 
point.  However  as  the  applied  displacement  u^^  gets  larger,  the 
error  becomes  larger  and  larger,  and  the  stress  and  strain  fields 
become  more  and  more  heterogeneous.  In  the  vicinity  of  the 
bifurcation  point  the  error  is  treated  as  an  artificial 
geometrical  imperfection  within  the  finite  element  mesh.  This 
artificial  imperfection  is  amplified  and  causes  heterogeneous 
stress  and  strain  solution.  This  geometrical  imperfection,  which 
comes  from  approximate  computation,  causes  the  same  effect  as  the 
introduction  of  an  eccentricity  in  the  problem  of  the  buckling  of 
an  Euler  beam.  A  difficulty  in  convergence  of  the  solution  is  also 


observed  when  the  applied  displacement  on  the  left  boundary  is  too 
large.  In  the  case  of  updated  Lagrangian  techniques  which  update 
nodes  position,  this  divergence  may  cause  the  annihilation  of  a 
few  elements  (negative  volume) ,  which  precedes  generally  a  fatal 
computational  crash. 

3.5  CHANGE  OF  EIGENVALUES  VERSUS  APPLIED  DISPLACEMENT. 

In  order  to  detect  the  bifurcation  point,  the  minimum 
eigenvalue  of  the  stiffness  matrix  is  calculated  at  each  step  of 
prescribed  displacement  u^.  The  minimum  eigenvalue  is  plotted  in 
Figs. 5  and  6.  versus  u^  for  three  different  meshes  and  for 
isotropic  and  anisotropic  materials.  As  shown  in  Figs.  5  and  6, 
the  minimum  eigenvalue  decreases  monotonously  and  becomes  zero  for 
some  value  of  the  applied  displacement  u^.  The  intersection  of 
A  .  with  the  u,  axis  corresponds  to  the  bifurcation  point,  where 
the  stiffness  matrix  becomes  singular.  As  shown  in  Figs. 5  and  6, 
the  mesh  refinement  has  a  very  important  influence  on  the 
detection  of  the  bifurcation  point.  The  analytical  value  of  the 
bifurcation  indicated  in  Table  3  is  overestimated  when  the  block 
is  not  enough  discretized.  However  25  elements  are  sufficient  to 
obtain  an  approximate  value  of  the  buckling  stress  of  the 
half-space  made  of  hypo-elastic  material. 

In  the  case  of  the  elastoplastic  model  of  Horii  and 

Nemat-Nasser,  1982,  the  minimum  eigenvalue  of  Fig. 7  remains 

positive  even  for  a  25  elements  model.  As  shown  in  Fig. 7,  first, 

A  .  starts  decreasing,  increases  when  u,  is  about  0.1,  and  then 
min  1 

decreases  again  for  larger  value  of  u, .  Since  A  .  does  not  become 
^  ^  1  min 

negative,  no  surface  instability  can  be  detected  for  this 
elastoplastic  material  during  the  wedge  test.  This  disagreement  of 
analytical  and  numerical  solutions  in  the  case  of  elastoplastic 
material  will  be  analyzed  in  a  later  section. 
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3.6  EIGENVALUES  AND  EIGENVECTORS  AT  THE  BIFURCATION  POINT. 


Once  the  bifurcation  point  has  been  detected  between  times  t 

n 

and  ^n+l'  iterations  may  be  performed  to  define  more  accurately 

the  bifurcation  time.  However  these  iterations  were  not  necessary 

in  our  analysis  since  the  prescribed  displacement  increment  was  as 

low  as  0.0025.  The  calculations  were  stopped  at  time  t^^^  when  the 

tangential  stiffness  matrix  was  found  to  have  a  negative 

eigenvalue.  The  eigenvalues  and  eigenvectors  of  the  stiffness 

matrix  were  calculated  at  the  onset  of  bifurcation  time  t  .  Tables 

n 

4  and  5  display  all  the  eigenvalues  in  the  case  of  4  elements  of 
isotropic  hypo-elastic  material  and  25  elements  of  anisotropic 
hypo-elastic  material.  The  eigenvalues,  which  are  sorted  according 
to  their  values,  are  listed  in  the  tables  4  and  5.  The 
corresponding  eigenvectors,  which  are  identified  by  the  mode 
number  of  the  eigenvalue,  are  plotted  in  Figs.  8  and  9.  The 
eigenvectors  are  plotted  with  respect  to  the  deformed  mesh,  which 
describes  the  position  of  the  body  at  the  onset  of  bifurcation. 
This  representation  of  eigenmodes  is  similar  to  the  one  used  for 
displaying  the  natural  modes  of  vibration  of  structural  systems. 
Figs. 8  and  9  show  the  direction  and  amplitude  of  the  increment  of 
nodal  displacement  which  are  an  alternate  solution  to  the  constant 
deformation  gradient  solution.  As  shown  in  Fig. 8a  to  8d  for  a  4 
elements  model,  the  modes  associated  with  the  lowest  eigenvalues 
are  surface  modes  for  which  the  amplitude  decays  with  depth.  The 
mode  of  Fig.8e,  which  corresponds  to  a  larger  eigenvalue,  is  a 
volume  mode.  Its  amplitude  does  not  decay  within  the  material.  In 
addition  to  the  eigenmodes  of  the  coarse  mesh  of  Figs. 8,  the 
eigenmodes  of  a  finer  25  elements  mesh  are  examined  in  Figs. 9.  The 
first  two  modes  of  Figs  9a  and  9b  are  symmetric  w.r.t  to  a  the 
central  vertical  axis.  The  amplitude  of  the  surface  wave  decays 
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Figure  9a.  First  eigenmode  for  a  25  efements  model  ot 
anisotropic  hy[)0-elastic  material. 


Anisotropic  hypo-olostlc  sotarlol  C"l.  Ct>0.  1  V>0. 3 


Figure  9d.  Fourth  eigenmode  for  a  25  elemen 
anisotropic  hypo-elastic  materia 


Fifth  eitjenmode  lor  a  25  elements  mode 
anisotropic  hypo-elastic  material. 
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ure  9g.  Seventh  eigenmode  tor  a  25  elcintmts  model  ol 
anisotropic  hy[)o-elast ic  material. 


Anisotropic  hypo-olostlc  aotorlol  C**t.  Gt**0.  1  V~0.  3 


Figure  9h.  Eitjth  eigeninode  for  a  2  5  elements  mcjde- 
anisotropic  hypo-elastic  mati-rial. 
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Figure  9i.  Ninth  eigenmode  for  a  elements  mode 
anisotropic  hypo-elastic  material. 


Anisotropic  hypo-alostlc  aotar 


(jure  9k.  Eleventh  eiqenmode  for  a  25  elements  iiujiie 
anisotropic  hypo-elastic  material. 


Anisotropic  hypo-«lastlc  ■oterlol  C~l.  Ct*0.  1  V>0.  3 


Fiqure  9n.  Fourteenth  eiqoniiiode  for  a  25  e-li'nient 
of  anisotrofjic  hypo-elastic  iiia  t  r  i  a  1  . 


Anisotropic  hypo-alostl 


F'iyure  9o.  I'ifteuiith  ei<jeninotle  for  a  2  5  (jliMuents  luudcJ 
anisotro[)ic  hypo-elastic  material. 


away  from  each  corner  in  the  and  directions,  which  was  not 
observed  in  the  case  of  the  half-space.  Following  these  two  first 
modes,  four  surface  modes  similar  to  the  ones  found  by  Biot,  1965, 
are  shown  in  Figs  9c  to  9f.  The  surface  wave  has  the  shortest  wave 
length  in  Fig. 9c.  This  minimum  wave  length,  which  is  associated 
with  the  lowest  eigenvalue,  is  directly  related  to  the  number  of 
nodes  on  the  free  surface.  Following  the  surface  modes,  four 
unexpected  modes,  which  are  not  mentioned  in  any  analytical  work, 
are  shown  in  Figs.  9g  to  9i.  These  unexpected  modes  are  associated 
with  edges  and  not  with  the  top  free  surface.  As  shown  in  Figs  9g 
to  9i  each  edge  of  the  block  has  a  similar  mode.  In  order  to 
designate  this  new  type  of  mode,  a  new  term  has  been  coined  -  edge 
mode.  Following  these  edge  modes,  the  volume  modes  are  shown  in 
Figs  9j  to  9p.  The  volume  modes  correspond  to  wavy  amplitude 
within  the  volume  of  the  block,  which  does  not  decrease  within  the 
material  volume.  Biot  (1965)  calls  these  volume  modes  internal 
buckling  modes.  The  correspondence  between  the  eigenvalues  of 
Table  5  and  the  types  of  mode  is  summarized  on  the  real  axis  of 
Fig. 10,  which  has  been  blown  up  between  0  and  1.  As  shown  in  Table 
5  and  Fig. 10,  all  the  eigenvalues  are  real.  The  minimum  and 
maximum  values  are  0.01397  and  14.28,  respectively.  In  a 
bifurcation  analysis,  all  these  modes  are  not  meaningful  and  only 
the  modes  with  the  eigenvalues  close  to  zero,  are  important. 
However,  as  it  is  shown  in  Fig. 6,  there  may  be  several  modes  with 
closely  spaced  eigenvalues,  which  may  define  several  possible 
bifurcated  branches  for  the  discretized  material.  This  multiple 
bifurcation  is  expected  in  the  case  of  surface  buckling.  Biot, 
1965,  found  that  there  is  an  infinite  number  of  surface  waves 
which  contribute  to  the  surface  buckling. 

In  the  case  of  the  anisotropic  hypo-elastic  material,  the  six 
first  surface  modes  may  contribute  to  the  instability.  According 
to  the  presence  and  the  relative  rank  of  surface  modes  with 
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respect  to  other  types  of  nodes,  it  may  be  concluded  that 
isotropic  and  anisotropic  hypo-elastic  materials  are  prone  to 
surface  instability  during  the  wedge  test.  However,  a  similar 
conclusion  cannot  be  reached  in  the  case  of  the  particular 
elastoplastic  model  which  was  selected  in  the  previous  section. 


Table  4.  Eigenvalues  of  the  stiffness  matrix  at  the  onset 
of  surface  instability  for  a  4  elements  models 
of  isotropic  hypo-elastic  material. 


Mode  No  Eigenvalues  Figure  No  Mode  No  Eigenvalues  Figure  No 


1 

-0.07849 

8a 

6 

3.315 

2 

0.2902 

8b 

7 

6.838 

3 

0.05912 

8c 

3 

9.320 

4 

0.8506 

8d 

9 

18.99 

5 

1.320 
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SECTION  4 


INFLUENCE  OF  MATERIAL  MODELLING  ON  SURFACE  INSTABILITY 


The  present  section  describes  the  influence  of  the  modelling 
of  material  behavior  on  the  surface  instability  of  the  wedge  test. 
The  following  constitutive  equations  have  been  used  to  simulate 
the  material  behavior: 

1.  fully  isotropic  hypo-elastic  materials,  with  two  material 
constants 

2.  transversely  isotropic  hypo-elastic  material,  with  three 
material  constants 

3 .  isotropic  elastoplastic  model  derived  from  the  flow  theory 
of  plasticity  with  five  material  constants 


The  parametric  study,  which  consisted  of  estimating  the 
influence  of  various  models  and  their  respective  material 
constants  on  the  buckling  load  of  the  wedge  test,  has  been  carried 
out  by  using  the  computer  program  BIF.  Before  presenting  the 
principles  of  the  program  BIF,  the  analysis  of  plane  strain 
surface  instability  is  summarized  in  the  form  introduced  by  Horii 
and  Nemat-Nasser  (1982).  Mathematical  derivations  can  be  found  in 
the  paper  by  Horii  and  Nemat-Nasser  (1982). 

4.1  PLANE-STRAIN  SURFACE  INSTABILITY  OF  HALF-SPACE. 


In  order  to  abbreviate  the  mathematical  notations,  it  is 
convenient  to  define  the  following  coefficients: 


'1111 


'2222 


^2 


(16a) 

(16b) 


5  9 


^1212 
"=1212  ^ 
^1212  ^ 


(16d) 

(16e) 

(16f) 

(16g) 

(16h) 


and  are  the  components  of  the  Kirchhoff  stress.  are 

the  coefficients  of  the  rate  type  constitutive  equations.  By  using 
the  coefficients  d^'s,  the  equilibrium  equations  of  incremental 


stress  in  plane  strain  condition  are: 


Vl,ll  ^  Vl,22  ^  (V^7)^2,12  =  ° 

'^5^2,11  ^  V2,22  ^  ° 


(17a) 


(17b) 


are  the  velocity  components.  The  subscript 


represents  spatial  derivative: 


k,  ij 


i » 3  » k— 1 f 2 


(17c) 


The  geometry  of  the  half-space  under  compressive  stress  is  shown 
in  Fig. 11.  It  is  supposed  to  be  infinite  in  the  vertical 
direction.  Note  that  the  coordinates  axis  x^-x^  have  been  replaced 
by  x-y  axis.  On  the  edges  x=  ±a,  a  uniform  stress  is  prescribed 
in  the  x  direction  without  shear  traction.  A  uniform  and  constant 
stress  is  applied  on  the  top  surface  (y=0) .  The  notion  of 

surface  instability  means  that  the  deformation  is  confined  close 
to  the  top  surface,  i.e.,  the  velocity  field  is  fading 

exponentially  with  x,  vanishing  at  infinite  distance  from  the  top 
surface.  Horii  and  Nemat-Nasser  (1982)  assumed  that  the  inception 
of  instable  deformation  is  defined  by  the  following  velocity 


field : 


-  iax  -aZy 
V,  =  A  e  .  e  ^ 


(18a) 


o  iaX  -aZy 
=  B  e  .e  ^ 


(18b) 


a  is  real  and  positive.  When  the  complex  variable  Z  has  a  real 
part  strictly  positive,  the  velocity  (v^/V^)  varies  as  cosine  or 
sine  in  the  x-direction  and  its  amplitude  decays  exponentially 
with  the  depth  y.  When  Z  is  purely  imaginary,  the  velocity 
varies  as  cosine  and  sine  in  the  x  and  y  directions  but  its 
amplitude  does  not  decay  with  depth.  The  former  case  corresponds 
to  surface  instability,  while  the  latter  case  corresponds  to 
volume  instability. 


Figure  11.  Half-space  under  compressive  stress. 


Since  the  velocity  field  must  decay  with  depth  for  surface 
instability,  the  velocity  field  becomes; 


The  coefficients  Z^,  Z^,  K^,  and  are  found  by  substituting  (19) 
into  (17).  Z^  and  Z^  are  the  complex  or  real  solutions  with 
positive  real  part  of  the  following  equation: 

aZ"^  -t-  bZ^  +  c  =  0  (19c) 

The  coefficients  of  the  quadratic  equation  (19c)  are: 


a  d2d2 

b  =  -  d3d5  -  d^d^  + 

c  -  d- dc 
1  0 


(19d) 

(19e) 

(19f) 


The  coefficients  and  are: 


-di  ^  Z,/d3 


^2  = 


-di  ^  Z2^d3 


(19g) 


The  imposition  of  the  boundary  conditions  of  Fig. 11  leads  to  two 
linear  equations  with  two  unknowns  and  A^: 


(idg  -  d^K^Z^  )A2  +  (  idg  -  d2K2Z2)A^  =  0 

(-dgZ^  +  iK^dJA^  +  (-d3Z2  +  iK^dJA^  -  0 


(20a) 

(20b) 


In  order  to  obtain  a  nontrivial  solution,  the  determinant  of  this 
system  must  vanish.  If  is  fixed  and  if  the  coefficients  d^  are 
assumed  to  be  function  of  r  only,  then  the  determinant  of  (20)  is 


a  function  of 


f(f^)  =  id2(Z^-Z2)(d5«lgd^-d2d^)  -Zj^Z2d2(d^^-d3dg;)  (21a) 

The  buckling  stress  of  the  half-space  can  be  obtained  by  solving 
the  following  nonlinear  equation  of  the  variable 

f(^l)  =  0  (21b) 

In  the  case  of  hypo-elastic  material,  all  coefficients  d^  are 
functions  of  r^.  However,  this  dependence  does  not  hold  for 
elastoplastic  models  due  to  the  influence  of  r^,  Horii  and 
Nemat-Nasser  (1982)  solved  this  problem  by  assuming  that  there  is 
no  plastic  flow  in  direction. 

4.2  COMPUTATION  OF  BUCKLING  LOAD  WITH  THE  PROGRAM  BIF. 

Based  upon  the  analytical  development  of  Horii  and 
Nemat-Nasser  (1982),  the  program  BIF  calculates  the  buckling 
stress  of  the  half  space  and  plots  the  corresponding  buckling 
eigenmode.  BIF  does  not  discretize  the  half  space  in  finite 
elements  but  solves  by  iteration  the  nonlinear  buckling  equation 
of  the  infinite  half  space  for  various  material  models.  That  is  to 
say  that  BIF  provides  analytical  results  which  can  be  compared  to 
our  finite  element  results.  In  order  to  solve  the  nonlinear 
equation  (21)  the  computer  program  increases  the  applied  stress 
in  a  segment  ^max^  which  is  specified  by  the  user.  Starting 

from  r  .  and  moving  by  increment  At  =  (r  ^  )/n  where  n  is  a 

specified  integer,  it  detects  the  stress  level  for  which  the  real 
or  imaginary  part  of  the  equation  (21)  changes  sign,  then  iterates 
by  using  Regula-Falsi  method  until  the  equation  (21)  is  satisfied. 
After  the  calculation  of  the  buckling  stress,  the  buckling 
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eigenmode  can  be  plotted  in  a  way  similar  to  finite  element 
post -processing . 

The  following  sections  summarize  for  each  model  the  influence 
of  the  material  constants  on  surface  instability. 

4.3  SURFACE  INSTABILITY  OF  ISOTROPIC  HYPO-ELASTIC  MODEL. 


Isotropic  hypo-elastic  material  have  two  material  constants: 
G  shear  modulus 

u  Poisson's  ratio 


If  the 

stress  on  the  free  surface  is  zero  t^=0, 

the  coefficients 

are: 

2G  +  A  —  T  ^ 

(22a) 

2G  +  A 

(22b) 

^  = 

G  -  ,^/2  = 

(22c) 

^5 

G  .  -  dj 

(22d) 

dg  =  A 

(22e) 

As  shown  in  Fig. 12,  the  buckling  stress  which  is 

normalized  by  the  shear  modulus  G,  increases  with  Poisson's  ratio 
and  depends  on  the  stress  applied  on  the  surface  of  the 
half-space.  As  it  is  to  be  expected,  compressive  stresses  acting 
on  the  half-space  surface  increase  the  buckling  load  while  tensile 
stresses  decrease  it.  Note  that  the  buckling  stress  is  larger  than 
the  shear  modulus  G,  which  is  unrealistic  for  most  Rock  Mechanics 
applications.  The  eigenmodes  or  buckling  modes,  which  describe  the 
velocity  field  at  the  inception  of  surface  instability,  have  been 
calculated  by  using  BIF.  The  mathematical  expression  of  these 
velocity  fields  can  be  found  in  Table. 6.  Three  buckling  modes  of 
the  half-space  are  plotted  in  Figs. 13a  to  13c  for  three  different 
values  of  the  wave  length  parameter  a  :  ir/2a,  3»r/2a,  and  5w/2a. 
The  eigenmodes  are  only  plotted  inside  a  square  window  of  depth  2a 
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which  is  centered  on  the  y-axis.  The  dashed  and  full  lines 
represent  the  positions  of  a  uniform  grid  before  and  after  the 
emergence  of  surface  instability,  respectively.  As  shown  in 
Figs, 13,  the  velocity  amplitude  decays  rapidly  with  the  depth  y. 
The  smaller  the  surface  wave  length  is,  the  faster  the  velocity 
amplitude  decays  with  depth.  Note  that  the  eigenmodes  of  the 
half-space  shown  in  Figs, 13  obey  almost  perfectly  the  boundary 
conditions  of  the  wedge  test  of  Fig, 2,  although  they  are  not 
required  to  do  so.  The  boundary  conditions  are  slightly  violated 
at  the  lower  boundary  of  Fig, 13a:  the  first  mode  of  the  half-space 
interferes  more  than  the  other  modes  with  the  boundary  conditions 
of  a  finite  block.  In  the  absence  of  an  analytical  solution  for 
the  surface  instability  of  block  of  finite  size,  this 

incompatibility  between  eigenmode  and  boundary  conditions  may 
explain  some  of  the  observations  made  during  the  finite  element 
analysis  of  the  wedge  test.  It  may  explain  why  the  surface  modes 
have  distinct  eigenvalues  and  why  the  lowest  eigenvalues  of  the 

stiffness  matrix  correspond  to  the  smallest  surface  wave  length. 

4.4  SURFACE  INSTABILITY  OF  ANISOTROPIC  HYPO-ELASTIC  MODEL. 

Anisotropic  hypo-elastic  have  three  material  constants: 

G,  shear  moduli  in  the  plane  of  isotropy  and  in  a  plane 

normal  to  the  plane  of  isotropy,  respectively 

1/  Poisson's  ratio 

As  shown  in  Fig. 14,  the  normalized  buckling  stress  decreases 
when  the  transverse  shear  modulus  G^  decreases.  The  case  G/G^  =  1 
corresponds  to  the  isotropic  case.  When  the  ratio  G/G^  becomes 
very  large,  which  corresponds  to  the  case  where  the  shear  modulus 

G  in  the  isotropy  plane  is  quite  larger  than  the  transverse  shear 

modulus  G^  in  a  plane  perpendicular  to  the  half-space  surface,  the 
buckling  stress  decreases  to  zero  but  remains  always  smaller  than 
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Figure  14.  Analytical  buckling  stress  of  anisotropic 

half  space  versus  ratio  of  shear  moduli  G/Gt  for 
Poisson  ratio  v  =  0  and  0.5. 


r.M: 


twice  G^.  The  buckling  modes  have  been  calculated  in  Table. 6  and 
plotted  in  Figs. 15  in  the  same  way  as  the  eigenmodes  of  the 
isotropic  case.  The  velocity  amplitude  of  Figs. 15a  and  15b  decays 
faster  with  depth  than  in  Figs. 13,  which  results  in  a  better 
compatibility  between  the  buckling  modes  of  the  half  space  and  the 
boundary  conditions  of  the  block  of  finite  size  of  Fig. 2.  Note 
that  the  modes  of  Fig. 15c  and  Fig. 13c  are  identical.  However,  the 
eigenmodes  of  the  half  space  of  Fig. 15  exhibit  a  larger  amplitude 
decay  with  depth  than  the  finite  element  eigenmodes  of  Figs  9c  to 
9f.  This  slight  difference  between  analytical  and  numerical  modes 
is  certainly  caused  by  the  spatial  discretization  used  in  the 
finite  element  model,  and  could  be  reduced  by  refining  the  finite 
element  mesh  in  the  vicinity  of  the  top  surface. 

4.5  SURFACE  INSTABILITY  FOR  ELASTOPLASTIC  MODEL. 

The  elastoplastic  model  used  in  the  analysis  is  an  isotropic 
elastoplastic  model  derived  from  the  flow  theory  of  plasticity.  It 
has  five  material  constants: 

G  elastic  shear  modulus 

1/  elastic  Poisson's  ratio 

n  friction  coefficient  i.e.  slope  of  yield  surface  in  the 

stress  space  (<7,r)  where  a  and  t  the  first  and  second 
deviatoric  stress  invariants,  respectively. 

0  dilatancy  parameter,  i.e.  slope  of  the  plastic  potential 
surface  in  the  stress  space  (a,r). 

H  plastic  modulus 

Further  detail  on  the  elastoplastic  models  may  be  found  in  the 
paper  of  Rudnicki  &  Rice  (1975),  or  Nemat-Nasser  &  Horii(1982). 
Nemat-Nasser  and  Horii  (1982)  applied  the  general  elastoplastic 
theory  to  describe  some  particular  experimental  data  obtained  on 
sandstones.  In  this  particular  model,  the  plastic  modulus  H  is  a 
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function  of  the  first  stress  invariant  a  and  the  effective  strain 
7.  The  friction  coefficient  /i  is  a  constant,  while  the  dilatancy 
angle  0  depends  on  the  stress  level  (see  Horii  and  Nemat-Nasser, 
1982,  for  further  detail). 

As  shown  in  Fig. 16,  the  buckling  stress  calculated  by  using 
elastic  plastic  model  depends  on  the  dilatancy  angle  0.  For  = 
0.39,  which  correspond  to  an  associative  flow  irule  since  n=0.29, 
the  normalized  buckling  stress  increases  with  the  normalized 

plastic  modulus  H/G.  In  the  same  way  as  the  eigenmodes  of 
hypo-elastic  models,  the  buckling  modes  of  elastoplastic  models 
have  been  calculated  in  Table. 6  and  plotted  in  Figs. 17  and  18  for 
two  values  of  the  normalized  plastic  modulus  H/G;  0.1  and  0.01. 
When  H/G=0.1,  the  velocity  amplitude  of  Figs. 17a  decays  with  depth 
but  certainly  not  as  quickly  as  in  Figs. 13  or  15.  The  eigenmodes 
of  Figs. 17  are  not  compatible  with  the  boundary  conditions  of  a 
block  of  finite  size,  when  H/G=0.01,  which  corresponds  to  more 
plastic  yielding  than  H/G=0.1,  the  eigenmodes  of  Figs. 18  are  found 
to  exhibit  a  very  slow  decay  of  amplitude  with  depth.  This  slow 
decay  with  depth  comes  from  the  small  real  part  of  the  root  Z^.  In 
the  case  H/G=0.01  the  eigenmodes  of  the  half-space  are  even  more 
incompatible  with  the  boundary  conditions  of  the  wedge  test.  The 
wave  length  parameter  a  must  be  very  large  to  compensate  for  the 
low  part  of  the  roots  and  Z^.  However,  in  contrast  to 

analytical  solutions,  finite  elements  are  unable  to  increase  the 
wave  length  parameter  a  in  order  to  compensate  for  the  slow 
amplitude  decay  resulting  from  the  small  real  parts  of  and  Z^, 
since  the  minimum  wave  length  of  the  discrete  system  is  controlled 
by  the  number  of  nodes  on  the  surface  as  it  was  pointed  out  in 
Fig. 9c.  The  combired  presence  of  a  small  amplitude  decay  and  a 
limited  wave  length  parameter  a  explains  why  the  surface 
instability  of  a  block  of  finite  size  is  difficult  to  obtain  in 
presence  of  large  plastic  yielding  within  a  discretized  system. 
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Figure  17c,  Third  eigenmode  of  the  half  space  for 
elastoplastic  model  (vi  =  3  =  0.39,  H/G=0. 
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Figure  18b.  Second  eigeninode  of  the  half  si^ace  for 
elastoplastic  model  ( ii  =  ti=^0 . 39  ,  11/Cj=0  .  1 
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The  influence  of  non-associative  flow  rule  has  also  been 
investigated  in  Figs. 19,  by  setting  the  dilatancy  parameter  ^ 
equal  to  zero.  In  this  case,  surface  instability  is  found  to  be 
only  possible  when  the  normalized  plastic  modulus  H/G  is  larger 
than  about  0.055.  For  value  smaller  than  0.055  the  roots  and 
of  equation  (20)  are  pure  imaginary,  which  implies  that  the 
amplitude  of  the  velocity  field  does  not  decay  with  depth.  The 
eigenmodes  corresponding  to  H/G=0.05  are  plotted  in  Figs. 19.  These 
eigenmodes,  which  show  no  sign  of  amplitude  decay  with  depth,  are 
incompatible  with  the  boundary  condition  of  the  wedge  test  and  of 
the  half-space.  From  the  results  of  Fig. 16,  it  may  be  concluded 
that  surface  instability  will  be  only  possible  when  stresses  reach 
a  value  in  the  order  of  the  elastic  shear  modulus,  which  is 
unlikely  to  take  place  due  to  plastic  yielding.  This  resistance  to 
surface  buckling  which  was  encountered  during  the  finite  element 
analysis  of  elastoplastic  material  subjected  to  the  wedge  test  was 
therefore  to  be  expected.  These  observations  are  in  agreement  with 
the  results  of  Hutchinson  and  Tvergaard  (1980)  but  in  disagreement 
with  the  results  of  Horii  and  Nemat-Nasser (1982) .  The  origins  of 
such  a  disagreement  need  to  be  investigated  in  more  detail  in 
future  research. 
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Elastic-plastic  madcl,  G=l.  V=0. 3.  Mu=0. 39.  Beta-0..  H=0, 05 


No . 1  -  Isotropic  hypo-elastic  (Figs. 13a  to  13c) 
G-l,  ~  I,  1/-0.3 

No . 2  -  Anisotropic  hypo-elastic  (Figs. 15a  to  15c) 
G-l,  -  0.1,  */  -0.3 

No. 3  -  Elastoplastic  model  (Figs. 17a  to  17c) 
G-l,  u  ~0. 3,  n  -  0  ~  0.39,  H-0.1 
No. 4  -  Elastoplastic  model  (Figs. 18a  to  18c) 
G-l,  1/-0.3,  ^  -  0  -  0.39,  H-0.01 
No . 5  -  Elastoplastic  model  (Figs. 19a  to  19c) 

G  -  1,  1/  -0.3,  fi  -  0.39.  0-0.,  H  -  0.05 
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SECTION  5 


INFLUENCE  OF  NONUNIFORM  STRESS  FIELD  ON  SURFACE  INSTABILITY 


The  fourth  task  of  the  project  was  to  examine  the  effect  of 
nonuniform  stress  distribution  on  the  value  of  the  buckling  load. 
This  influence  has  been  investigated  by  restraining  the  vertical 
displacement  of  the  nodes  located  on  the  left  boundary  of  Fig. 2. 
The  nodal  displacement  are  prescribed  in  both  horizontal  and 
vertical  direction.  This  prescribed  displacement  boundary 
condition  corresponds  to  the  case  of  no  slippage  between  the 
material  and  the  left  loading  plate.  Besides  the  vertical 
displacement  conditions  on  the  left  boundary,  the  geometry  and 
boundary  conditions  of  the  model  are  similar  to  the  ones  shown  in 
Fig. 2.  The  deformed  mesh  at  the  onset  of  surface  instability  is 
shown  in  Fig. 20.  As  shown  in  Fig. 21,  the  stress-displacement 
response  is  almost  identical  to  the  frictionless  case;  however  the 
introduction  of  friction  tends  to  make  the  material  response 
slightly  stiffer  than  in  the  absence  of  friction.  Figs. 22a,  b,  and 
c  represent  the  contour  of  the  various  components  of  stress 

and  which  characterize  the  degree  of  nonuniformity  of  the 
stress  field.  The  horizontal  stress  varies  between  -0.270  and 
-0.241,  the  stress  between  -0.058  and  +0.0036,  and  the  shear 
stress  £7^2  between  -0.00182  and  0.00629.  All  stresses  can  be 
considered  dimensionless  since  the  shear  modulus  G  has  been 
selected  equal  to  1.  The  detection  of  the  bifurcation  can  be 
performed  in  the  same  way  as  in  the  absence  of  friction.  The 
variation  of  the  minimum  eigenvalues  of  the  stiffness  matrix 
versus  the  applied  displacement  is  shown  in  Fig. 23.  This  figure 
also  indicates  that  the  minimum  eigenvalue  of  the  tangential 
stiffness  matrix  with  friction  is  always  larger  that  the  one 
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on  the  left  bcjundary. 


acement  responses  with  and  without 
anisotropic  hyijo-eiast  ic  material 


Figure  22b.  Contours  ot  stress,  in  presence  o 

for  anisotropic  hyijo-olastic  material. 


without  friction.  Therefore  introducing  friction  on  the  boundary 
increases  the  buckling  load.  The  eigenvectors  associated  with  the 
eigenvalues  of  Table  7  are  shown  in  Figs. 24a  to  24i.  The 
eigenmodes  are  similar  to  the  ones  obtained  in  the  absence  of 
friction,  with  the  difference  that  the  amplitude  on  the  left 
boundary  are  zero.  The  mode  associated  with  the  lowest  eigenvalue 
is  quite  similar  to  the  second  mode  which  was  obtained  for  the 
frictionless  conditions  of  Figs. 9.  The  similitude  pertains  for  the 
surface  modes  of  Figs.  24b  to  24e  and  for  the  edge  modes  of  Figs. 
24f  to  24h.  Two  edge  modes  of  Figs. 9  are  missing  due  to  the 
kinematic  constraints  imposed  at  the  left  boundary  of  Fig. 2.  The 
buckling  load  increases  as  friction  is  introduced  since  the 
kinematic  constraints  prevent  the  buckling  modes  from  being 
activated. 
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Table  7 


Eigenvalues  of  the  stiffness  matrix  at  the  onset 
of  surface  instability  for  a  25  elements  model 
of  anisotropic  hypo-elastic  material  in  presence 
of  friction. 


Mode  No 


Eigenvalue  Type 


Mode  No  Eigenvalue 


1 

0.01455  surface 

26 

3.008 

2 

0.04254 

27 

3 . 390 

3 

0.08934 

28 

3.878 

4 

0.1359 

29 

3 .958 

5 

0 .2274 

30 

4.415 

6 

0.4679  edge 

31 

4.815 

7 

0.6216 

32 

4.872 

8 

0.6433  volume 

33 

4.872 

9 

0.7814 

34 

5.203 

10 

1.030 

35 

5.465 

11 

1.039 

36 

5.715 

12 

1.171 

37 

5.746 

13 

1.268 

38 

6.923 

14 

1.305 

39 

7.100 

15 

1.574 

40 

7.153 

16 

1.584 

41 

7.546 

17 

1.672 

42 

8.392 

18 

1.959 

43 

9.170 

19 

2.013 

44 

9.239 

20 

2 . 088 

45 

9.760 

21 

2 . 168 

46 

10.54 

22 

2.449 

47 

10.84 

23 

2 . 682 

48 

12.49 

24 

25 

2.768 

2 . 798 

49 

14.19 

SECTION  6 


CONCLUSIONS  AND  RECOMMENDATIONS 


6.1  CONCLUSIONS. 

A  general  numerical  approach  has  been  developed  to  assess  the 
phenomenon  of  roc)cbursting ,  based  upon  the  assumption  that 
rockbursting  is  a  surface  instability  of  bulky  material  masses. 
The  general  approach  has  been  calibrated  in  the  case  of  a  simple 
boundary  value  problem:  the  wedge  test. 

The  analysis  has  been  carried  out  in  four  steps: 

1.  Using  nonlinear  finite  element  methods  similar  to  the  ones 
employed  in  structural  dynamics,  finite  element  computer 
program  modules  have  been  developed  for  analyzing  surface 
instability  of  material  masses. 

2.  The  modules  have  been  used  to  analyze  a  simple  boundary  value 
problem:  the  wedge  test.  The  numerical  solutions  have  been 
compared  with  analytical  solutions  in  the  case  of  isotropic  and 
anisotropic  hypo-elastic  materials. 

3.  Three  constitutive  equations  formulated  in  terms  of  large 

deformation  have  been  used  to  describe  the  material  behavior: 
isotropic  and  anisotropic  hypo-elastic,  and  elastoplastic 

models.  The  influence  of  their  respective  material  constants  on 
the  surface  instability  during  the  wedge  test  has  been 
evaluated. 

4.  The  influence  of  the  friction  on  surface  instability  has  also 
been  estimated  in  the  case  of  the  wedge  test. 

Surface  instability  has  been  treated  as  a  bifurcation 

phenomenon,  which  emerges  from  the  loss  of  uniqueness  of  the 
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partial  differential  equation  of  incremental  stress  equilibrium. 
In  the  context  of  finite  elements,  the  bifurcation  point  is 
detected  when  the  tangential  stiffness  matrix  becomes  singular. 
Accurate  finite  element  methods  have  been  required  in  order  to 
detect  the  bifurcation  and  to  avoid  numerical  divergence.  In  the 
vicinity  of  a  bifurcation  point,  finite  element  methods  have  been 
found  to  consider  the  computational  errors,  which  are  inherent  due 
to  any  iterative  techniques,  as  a  geometrical  imperfection  of  the 
finite  element  mesh.  In  some  circumstances,  these  artificial 
imperfections  were  found  to  grow  rapidly  and  to  transform  the 
uniform  stress-strain  state  into  highly  heterogeneous  stress  and 
strain  field  leading  to  numerical  divergence.  These  artificial 
geometrical  imperfections  have  the  same  effect  as  the  eccentricity 
in  the  buckling  problem  of  Euler  beams.  From  another  point  of 
view,  these  artificial  imperfections  render  difficult  the  approach 
of  bifurcation  point  by  transporting  the  numerical  solution  on  a 
bifurcated  branch  of  solution  far  away  from  the  analytical 
bifurcation  point. 

Isotropic  and  anisotropic  hypo-elastic  models  were  found  to  be 
simple  to  use  in  estimating  the  buckling  load  and  eigenmode  of  the 
wedge  test.  For  such  hypo-elastic  materials,  numerical  and 
analytical  solutions  were  in  excellent  agreement,  which 
demonstrated  that  our  numerical  methods  and  implementation  were 
appropriate  to  describe  surface  instability.  The  elastoplastic 
models  of  the  flow  theory  of  plasticity  which  are  used  in  this 
analysis  are  found  to  be  resistant  to  surface  instability  during 
the  finite  element  analysis  of  the  wedge  test.  A  more  detailed 
analytical  investigation  corroborated  these  finite  element 
observations.  This  analysis  revealed  that  a  slice  of  the  half 
space  made  of  the  same  elastoplastic  models  is  prone  to  volume 
instability  for  low  stress  level  and  to  surface  instability  for 
high  stress  level.  In  other  words,  this  analysis  pointed  out  that 


surface  instability  may  be  impossible  for  some  combinations  of 
material  parameters  or  loading  conditions  since  plastic  yielding 
prevents  stress  from  reaching  values  of  the  same  order  as  elastic 
shear  modulus.  These  observations  on  the  influence  of  material 
modeling  on  surface  instability  are  in  agreement  with  Hutchinson 
and  Tvergaard  conclusions  (1980)  but  in  disagreement  with 
Nemat-Nasser  and  Horii  results  (1982). 
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6.2  RECOMMENDATIONS. 


In  order  to  clarify  the  problems  encountered  in  the  wedge  test 
of  elastoplastic  model,  it  is  recommended  to  extend  the  methods 
used  in  the  program  BIF  from  the  half-space  to  block  of  finite 
length.  This  extension  will  provide  a  valuable  approach  for  the 
investigation  of  surface  instability  in  the  case  of  block  of 
finite  size  made  of  elastoplastic  material. 

The  numerical  methods  which  have  been  developed  in  this  work 
have  only  been  applied  to  the  simple  boundary  value  problem  of  the 
wedge  test.  According  to  the  good  agreement  between  numerical  and 
analytical  results  in  the  case  of  the  wedge  test,  it  is 
recommended  to  extend  the  method  to  the  more  complex  boundary 
problems  encountered  in  real  excavations.  A  direct  application  of 
this  generalization  is  to  estimate  the  risk  of  rockbursting  of 
real  excavations. 

The  present  work  analyzed  only  the  emergence  of  surface 
instability  but  did  not  examine  how  the  phenomenon  evolve 
afterwards.  In  order  to  extend  the  solution  in  the 
post-bifurcation  range,  it  is  recommended  to  use  the  computational 
techniques  of  structural  mechanics  known  as  continuation 
techniques  or  imperfection  approach. 

Finally  it  is  suggested  to  perform  careful  experiments  on 
material  prone  to  surface  instability  and  to  develop  constitutive 
models  which  are  more  realistic  than  the  ones  used  in  this  work. 
Extreme  care  must  be  exercised  in  the  development  of  such  material 
models.  The  anisotropic  hypo-elastic  model  used  in  this  analysis 
has  clearly  indicated  that  surface  instability  may  be  dependent  on 
material  behavior  which  are  not  apparent  in  the  stress-strain 
response  prior  to  the  buckling  phenomenon.  It  is  recommended  to 
develop  models  based  upon  vertex  plasticity  or  the  deformation 
theory  of  plasticity. 
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APPENDIX  A 


I 

BIOT'S  APPROACH  TO  SURFACE  INSTABILITY 


The  following  section  presents  Biot's  approach  to  investigate 
the  problem  of  the  surface  instability  of  an  elastic  half-space 
subjected  to  finite  deformation.  It  is  a  simple  illustration  for 
the  understanding  of  the  phenomenon  of  surface  instability  from 
the  point  of  view  of  continuum  mechanics.  The  problem  was 
initially  formulated  by  Biot  in  his  book  "Mechanics  of  incremental 
deformation,"  John  Wiley  &  Sons,  Inc.  1965,  page  159.  This  pioneer 
work  defines  the  problem  of  surface  instability  in  a  very  simple 
but  rigorous  way. 

Biot  describes  the  material  behavior  with  incompressible 
finite  deformation  elasticity: 

a  =  MqG  -  p.S  (la) 

a  is  Cauchy  stress,  Mq  the  only  material  constant,  p  the  mean 
pressure,  S  is  Kroneker's  symbol  with  components 

=  1,  if  i=j ;  and  5^,^  =  0,  if  i^-j  (lb) 

G  is  related  to  the  deformation  gradient  F: 

G  =  P.F^  F.  .  =  (2) 

aXj 

where  x  and  X  are  the  initial  and  present  positions  of  a  particle, 
respectively.  The  constitutive  equation  (1)  can  be  written  in 
incremental  form: 


a  —  (D*G'^G*D)  —  p*6 


D  is  the  rate  of  deformation: 


dU.  dU . 

D.  .  =  (  —  +  — L  )/2 
^  axj  dx^ 


where  u  is  the  velocity  at  location  x 


"ij  =  ^ijkl  °kl 


a,  Jauman's  rate  and,  a  material  time  derivative  of  Cauchy  stress 
a  are  related  (see  page  32  of  Biot's  book,  1965) 


-  CT  .  W  +  W.  £7 


+  a  • .  W.,  +  O',  W., 

Ik  3k  3k  Ik 


W  is  the  spin  tensor 


au .  au. 

w..  =  (  ^  -  -L  )/2 

ax^  ax^ 


The  incremental  constitutive  matrix  in  the  relation  (3)  is 


^ijkl  “  ^  ^ij^kl  ^  ^®ik^lj  ®il^jk  ®jl^ik  ^  ®jk^il^/^ 

s.  .  is  the  deviator  stress 
11 

®ij  =  ^ij  ■  W^^ij 


'J'  'V*  "w"  •"  ^  I  ^  **-  •’*-  -‘k  ^•*  •**  ^  -  **  •  ”*• 


tf- 


A  is  an  arbitrary  large  number,  which  accounts  for 
incompressibility. 

The  problem  formulated  by  Biot  is  shown  in  Fig. 25. 


X  =x 

(A\) 


A 

2 

Figure  25.  Half-space  under  compressive  stress. 

The  trivial  solution  of  the  problem  of  Fig.  25  is  a  unifoirm 
displacement  field,  with  homogeneous  strain.  According  to  the 
material  incompressibility  and  plane  strain  condition,  the 
stretches  in  the  1,  2  and  3  directions,  A^,  A^,  and  a^  are  such 
as : 


Using  the  constitutive  equation 


(9) 

(1),  the  stresses  are 

(10) 

(11) 
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By  substracting  (10)  from  (11),  the  applied  stress  P  is; 

Using  the  incremental  constitutive  relation  (3) ,  the  increment  of 
convected  stress,  which  is  Jauman's  rate  of  Cauchy  stress,  is 


.11  -  ^  = 

2. 

(14) 

^22  -  ^  = 

2m  0^2 

(15) 

A 

"12  = 

2m  D^2 

(16) 

Since  the  material  is  incompressible  °  found  by 

adding  (14)  and  (15): 


(17) 


The  modulus  n  is 


^  +  ^2  )/2 


(18) 


Biot  derives  the  equilibrium  conditions  for  the  incremental  stress 
field: 


aw 


-  p 


12 


=  0 


(19) 


ax 

ay 

ay 

da 

¥  — — — 

aw, 

-  p  — 

12 


=  0 


(20) 


ax 


ay 


ax 
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The  constraint  of  incompressibility  is  satisfied  by  introducing  a 
stream  function  «j(x,y): 


dip 

dip 

— 

u  =  — 

(21) 

ay 

^  ax 

After  some  calculations,  Biot  found  equations  (19)  and  (20) 
become : 


5  <p  d  ((> 

(;i-P/2)  —  +  (/i+P/2)  -y 
3x  ay 


0 


(22) 


The  usual  solution  of  (22)  is  the  uniform  displacement  field  with 
homogeneous  deformation.  However  there  is  a  solution  which  is 
sinusoidal  along  the  x  direction  and  vanishes  at  infinite  depth 
(y=-oc)  : 

V  =  -^(C^e^^  +  C2e^-^^)sin  lx  (23) 

and  are  two  constants,  1  is  inversely  proportional  to 

length,  and  k  is 


k  -  (24) 

By  enforcing  the  fact  that  the  surface  y=0  is  free  of  stress,  Biot 
finds  that  the  solution  (23)  exists  when: 

-  2  =  0  with  f  =  (l-k^)/(l+)c^)  (25) 

which  has  one  real  solution  f  =  0.839  (k=0.295).  The  extension 
ratio  corresponding  to  surface  instability  for  this  case  is 

=  0.544  (26) 
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and  the  applied  compressive  stress  is 


P  =  3 . 08  Mr 


APPENDIX  B 


ANALYTICAL  SOLUTION  OF  THE  RESPONSE  OF  AN  HYPO-ELASTIC  SOLID 
SUBJECTED  TO  FINITE  PLANE  STRAIN  COMPRESSION 


Consider  the  plane  strain  problem  of  Fig. 26, 


liiiliiiiii  a. 


Fig. 26.  Block  of  finite  size  subjected  to  the  wedge  test, 


The  constraints  imposed  by  the  plane  strain  loading  on  the  rate  of 
deformation  are: 

^33  ^  °12  "  °23  "  °31  ^  ^ 

=  0  i,j=l,2,3  (lb) 

The  constitutive  equation  is  expressed  in  term  of  the  Jauman  rate 
of  Cauchy  stress: 

a  =  C.D  -  a  trace(D)  (Ic) 


lVS  ^ 


where  a  is  Cauchy  stress,  and  C  the  constitutive  matrix  of 
transversely  anisotropic  elasticity,  the  plane  of  isotropy  being 
normal  to  the  direction.  According  to  (lb)  Jauman's  rate  is 

identical  to  the  material  time  derivative 


a 


a 


(Id) 


The  general  constitutive  (la)  becomes; 


11 


2G 


l-2u 


((!■")  Dll  ^  "D22)  -  ^ii(Dii-^D22) 


(2a) 


22 


2G 


l-2u 


("□ii  +  (1-^)D22)  -  ^22(Dii+D22) 


(2b) 


33 


2G 


l~2u 


(i/(Dj_^+D22)  )  " 


(2C) 


^^12  "  ^t  Di2  -  D 


(2d) 


If  the  stress  on  the  free  surface  X2=0  remains  equal  to  zero: 

O 

£722  =  0  and  <722  =  0  (3a) 

then  the  relation  (2b)  reduces  to: 


D22 


1-1 


11 


and  the  relation  (2a)  becomes 


(3b) 
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^  /  / 


The  relation  (4)  can  be  integrated  with  respect  to  time  t  assuming 
that  =  0  and  =  0  at  time  t  =0 


D  dt  - - Ln^  ^  ) 

0  l-2u  2G 


But  the  integral  of  is  related  to  the  displacement  in  the 

x^-direction  u^  which  is  zero  at  time  t=0 


dt  =  Ln(l+u^) 


The  material  response  is: 


l-2u 


l-2i/ 


*^33  ‘'ll 


[1  -  (1+u^) 


u^  =  -  1  +  (1+U^) 


Note  that  the  material  response  is  independent  of  the  modulus 
since  there  is  no  shear.  The  finite  strain  solution  (5)  is  equal 
to  the  snail  strain  solution  if  u.  is  small  compared  to  1 


u^  =  -^u^/(l-i/) 


response  given  by  equation  (5)  is  plotted  in  Figs- 3  and  27  fo 
following  values  of  the  shear  modulus  and  Poisson's  ratio 
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